Aberrant expression detection in ATG7 patients
David Zhang
Aim: Detect pathways/genes that are differentially expressed in response to disease.
1 Background
1.1 Patient metadata
The 5 patient samples of interest are clinically diagnosed with mitochondrial disease and genetically diagnosed with mutations in the gene ATG7. Details below:
load(here::here("results/aberrant_expression/OUTRIDER/50_controls/gene_count_ods.rda"))
gene_count_ods <- gene_count_ods[, colData(gene_count_ods)[["samp_id_tidy"]] != "S2557"]
which_ATG7 <- which(colData(gene_count_ods)[["gene_name"]] == "ATG7")
which_ATG7 <- colData(gene_count_ods)[["samp_id_tidy"]][which_ATG7] %>%
as.character()
ATG7_patient_metadata <- colData(gene_count_ods)[which_ATG7,] %>%
as.tibble() %>%
dplyr::select(samp_id_tidy, gene_name, mutation_type, mutation_hgvs)
ATG7_patient_metadata$mutation_type <- ATG7_patient_metadata$mutation_type %>% lapply(FUN = str_c, collapse = ";") %>%
unlist()
ATG7_patient_metadata$mutation_hgvs <- ATG7_patient_metadata$mutation_hgvs %>% lapply(FUN = str_c, collapse = ";") %>% unlist()
ATG7_patient_metadata <- ATG7_patient_metadata %>%
mutate(family = c(2, 2, 3, 4, 1)) %>%
arrange(family)
ATG7_patient_metadata %>%
data.table() ATG7_patient_metadata %>%
write_delim(here::here("results/tables_figures/ATG7_patient_metadata.csv"), delim = ",")1.2 Publications
- RNA-SEQC - https://academic.oup.com/bioinformatics/article/28/11/1530/267467
- OUTRIDER - https://www.sciencedirect.com/science/article/pii/S0002929718304014?via%3Dihub
- DESeq2 - https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8
- megadepth - https://github.com/ChristopherWilks/megadepth
- dasper - https://github.com/dzhang32/dasper
- gprofiler2 - https://academic.oup.com/nar/article/47/W1/W191/5486750
1.3 OUTRIDER
- OUTRIDER is a bioinformatics software that aims at detecting aberrantly expressed genes from RNA-seq data. It is specifically designed for a rare disease context (1 vs all) by correcting for co-variations that are not known a-priori using an autoencoder to to help distinguish outliers more readily. OUTRIDER then applies a negative binomial test to find outlier expressed genes.
1.4 DESeq2
- DESeq2 aims at finding differentially expressed genes by modelling reads using a negative binomial distribution (similar to OUTRIDER), then using a binomial test to ascertain significance.
2 Method
Patient and control samples
RNA-sequencing was performed on a total of 58 individuals. 5 of these were ATG7 patients detailed in table X. In particular, samples were obtained from patients M1856.17, M0920.18, M0921.18, M1111.19 and M1716.19 from families 1, 2, 2, 3 and 4 respectively. The remaining 53 samples were used as the control samples consisting of a cohort of 46 patients diagnosed with mitochondrial disease, 7 patients with congenital muscular dystrophyand 2 unaffected individuals.
Culturing fibroblasts
Fibroblast cell lines were obtained from 58 individuals were obtained through collaboration with Dr. Haiyan Zhou and the Lily Consortium consisting of Dr. Charu Deshpande, Dr. Ines Barbosa, Prof. Joanna Poulton, Prof. Michael Simpson, Prof. Robert McFarland, Dr. Robert Pitceathly and Prof. Robert Taylor. DMEM supplemented with 10% Fetal Bovine Serum and 0.05 g/ml Uridine was used for culturing cells. Harvesting cells was performed by detaching cells using TrypLE Enzyme, washing cells once with washed with DPBS before storage at -80°C.
RNA-sequencing
RNA from fibroblast pellets was either extracted and sequenced at Eurofins Genomics or extracted through Bioxpedia, then sequenced at UCL Genomics. RNA integrity numbers (RIN) were measured using Agilent Technologies 2100 Bioanalyzer or Agilent 4200 Tapestation by Eurofins/Bioxpedia respectively. All RIN values exceeded 8.0. RNA library preparation was performed using the INVIEW Transcriptome Discover protocol from Eurofins. Specifically, the cDNA library was prepared using a random primed, strand specific, poly-A selection protocol. Paired-end sequencing was performed on illumina NovaSeq 6000 Sequencing System machines to an coverage of ~100 million pair-end reads per sample. Reads of 150-bp length were used to increase the proportion of reads with a gapped alignment, which represent splicing events.
RNA-sequencing alignment and quality control of patient samples
Pre-alignment quality control including adapter trimming, read filtering and base correction was performed using fastp, an all-in-one FASTQ preprocessor (v0.20.0) (CITATION - https://academic.oup.com/bioinformatics/article/34/17/i884/5093234). Reads were aligned using STAR 2-pass (v2.7.0) to the hg38 build of the reference human genome (hg38) using gene annotation from Ensembl v97 (CITATION - https://pubmed.ncbi.nlm.nih.gov/23104886/). 1st pass alignment was used to discover novel junctions, which were used as input to the 2nd pass to improve the sensitivity of junction detection. Via the option --outFilterMultimapNmax 1 reads were required to be uniquely mapped to only a single position in the genome. The minimum required overhang length of an annotated and unannotated junction was required to be 8 and 3 base pairs respectively. The output BAM files underwent post-alignment QC using RSeQC (CITATION - https://academic.oup.com/bioinformatics/article/28/11/1530/267467) with all samples passing quality control after manual assessment.
Detecting abberantly expressed genes and pathways
Gene count matrices were generated using RNA-SeQC (v2.3.4) with Ensembl v97 reference annotation matching the protocol used in https://github.com/broadinstitute/gtex-pipeline (CITATION - https://pubmed.ncbi.nlm.nih.gov/22539670/). OUTRIDER was used to find outlier expression of genes in each ATG7 patient, hence an all vs 1 experiment design was used, whereby each ATG7 patient were compared to the remaining 53 controls (CITATION - https://www.sciencedirect.com/science/article/pii/S0002929718304014?via%3Dihub). OUTRIDER was run with default setting and the autoencoder was fit using the optimal encoding dimension found to the be XX. P-values were corrected for multiple testing using the Benjamin Hochberg method and 0.05 was used as the signifance threshold. DESeq2 was applied to detect common genes or pathways that were dysregulated across all ATG7 patients (CITATION - https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8). Thus, all 5 ATG7 were grouped together as the case cohort and compared to the 53 control samples. The Benjamin Hochberg method was used for multiple test correction and genes with a p-value less than 0.05 and an absolute log-fold change greater than 1.5 were considered differentially expressed. gprofiler2 was then run to with differentially expressed genes as input, ranked by p-value to find dysregulated pathways (CITATION - https://academic.oup.com/nar/article/47/W1/W191/5486750).
Detecting abberantly expressed splicing events
dasper (https://github.com/dzhang32/dasper) was used to detect aberrant splicing events in each of the 5 ATG7 patients. Foremost, BAM files were converted into the BigWig format using megadepth (CITATION - https://github.com/ChristopherWilks/megadepth). Junctions outputted by STAR and BigWig files were inputted into the dasper pipeline. Foremost, dasper loads, filters junction data. Here, we required junctions to have at least 5 counts in at least 1 sample, a length between 20-1,000,000 base pairs and to not overlap any ENCODE blacklist regions (CITATION - https://www.nature.com/articles/s41598-019-45839-z). Junctions were then annotated using Ensembl v97 reference annotation and classified in the categories "annotated", "novel_acceptor", "novel_donor", "novel_combo", "novel_exon_skip", "ambig_gene" and "unannotated". Junctions were then clustered by their shared acceptor or donor sites and their counts were normalised. Coverage data associated with each junction was loaded and normalised for each sample. Junction and coverage normalised counts were scored using the z-score approach. Resulting z-scores for each patient were placed into a isolation forest model to rank splicing by how abberrant they looked in relation to controls. The rank for the most abnormal splicing event in ATG7 was then extracted. Sashimi plots describing the splicing across the ATG7 patients were plotted for each of the 5 patients.
2.1 OUTRIDER methods
2.1.1 Filter by expression
- There are 12339 genes remaining after filtering by expression.
knitr::include_graphics(here::here("results/aberrant_expression/OUTRIDER/50_controls/FPKM.png"))
knitr::include_graphics(here::here("results/aberrant_expression/OUTRIDER/50_controls/expressed_genes.png"))
2.1.2 Pre-correction gene count heatmaps
knitr::include_graphics(here::here("results//aberrant_expression/OUTRIDER/50_controls/count_cor_pre_correction.png"))
knitr::include_graphics(here::here("results//aberrant_expression/OUTRIDER/50_controls/count_gene_sample_pre_correction.png"))
2.1.3 Optimising the autoencoder correction dimension
knitr::include_graphics(here::here("results//aberrant_expression/OUTRIDER/50_controls/encoding_dim.png"))
2.1.4 Post-correction gene count heatmaps
Here, I have included one example (M0920.18). In reality, there's plots like these for each ATG7 sample as I have re-run each ATG7 with the other ATG7 samples removed.
knitr::include_graphics(here::here("results//aberrant_expression/OUTRIDER/50_controls/count_cor_post_correction_M0920.18.png"))
knitr::include_graphics(here::here("results//aberrant_expression/OUTRIDER/50_controls/count_gene_sample_post_correction_M0920.18.png"))
3 Results
3.1 Obtain gene counts
3.1.1 Generate genes-only version of Ensembl v97 gtf
cd /tools
git clone https://github.com/broadinstitute/gtex-pipeline
chmod -R 775 gtex-pipeline
# generate the genes gtf
/tools/gtex-pipeline/gene_model/collapse_annotation.py \
/data/references/ensembl/gtf_gff3/v97/Homo_sapiens.GRCh38.97.gtf \
/data/references/ensembl/gtf_gff3/v97/Homo_sapiens.GRCh38.97.genes.gtf3.1.2 RNA-SEQC
source(here::here("scripts/aberrant_expression/get_gene_count_RSE.R"))3.2 Which pathways/genes are differentially expressed in response to disease?
3.2.1 ATG7 expression is significantly lower in family 1
We compared the expression of ATG7 measured as fragments per kilobase of exon model per million reads mapped (FPKM) against the remaining patient samples. Family 1 was observed to have lower ATG7 expression compared to the other patients. All other families looked to have similar ATG7 expression to the remaining patient samples.
plot_fpkm <- function(gene_id, gene_count_ods, ATG7_patient_metadata){
gene_count_ods_filtered <- gene_count_ods[names(gene_count_ods) == gene_id, ]
expr_to_plot <-
assays(gene_count_ods_filtered)[["fpkm"]] %>%
as.data.frame() %>%
gather(key = "samp_id", value = "fpkm") %>%
mutate(ATG7_patient = samp_id %in% ATG7_patient_metadata[["samp_id_tidy"]])
families <- tibble(samp_id = c("M0920.18", "M0921.18", "M1111.19", "M1716.19", "M1856.17", "S2557"),
family = c(2, 2, 3, 4, 1, 1))
expr_plot <- expr_to_plot %>%
left_join(families) %>%
mutate(family = ifelse(is.na(family), "non-ATG7 sample", family)) %>%
ggplot(aes(x = ATG7_patient, y = fpkm)) +
geom_boxplot() +
geom_jitter(aes(colour = family)) +
scale_x_discrete(name = "ATG7 Patient") +
scale_y_continuous(name = "FPKM") +
scale_colour_manual(values = c(ggpubr::get_palette("npg", 4), "black")) +
theme_bw()
return(expr_plot)
}
ATG7_expr_plot <- plot_fpkm(gene_id = "ENSG00000197548",
gene_count_ods,
ATG7_patient_metadata)
ATG7_expr_plot
ggsave(ATG7_expr_plot,
filename = "ATG7_expr_plot.png",
path = here::here("results/tables_figures"),
width = 8, height = 6, dpi = 600)3.2.2 ATG7 is detected an expression outlier in family 1
Consistent with the expression plot above, OUTRIDER detected ATG7 as an expression outlier in family 1 with a negative z-score of -4.59.
load(here::here("results/aberrant_expression/OUTRIDER/50_controls/ATG7_res_all.rda"))
ATG7_res_ATG7_only <- ATG7_res_all %>%
filter(sampleID %in% ATG7_patient_metadata[["samp_id_tidy"]]) %>%
mutate(padjust = p.adjust(pValue, method = "fdr")) %>%
dplyr::select(Description, everything())
ATG7_res_ATG7_only %>%
filter(Description == "ATG7")3.2.3 OUTRIDER detects 9 genes that are expression outliers common to at least 2 of the ATG7 patients
OUTRIDER detected 264 unique genes as expression outliers across all samples (Table ATG7_OUTRIDER_res_signif.csv). In order to find genes that were commonly disrupted across ATG7 patients we filtered results by genes that were outliers in at least two ATG7 patient samples. The table below details the 9 genes left after this filtering. Each row represents a patient and each column represents a gene that was detected as a outlier in at least 2 of the ATG7 patients. Values in each cell represent the “z-score;pvalue” for each patient/gene respectively. Green indicates a positive z-score and red a negative z-score with respect to controls. Bold font is used for those that are significant outliers (p-value < 0.05 after FDR correction).
VPS41, HMG20A and TBC1D3L were specifically down-regulated in family 2. TBC1D3L is a GTPase activating protein for RAB5, which looks to have a role in macroautophagy (https://jcs.biologists.org/content/121/10/e1002). VPS41 has a role in vesicular transport and is suggested to be neuroprotective in AD/PD acting via an autophagy mechanism (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6276837/).
ATG7_res_ATG7_only_signif <- ATG7_res_ATG7_only %>%
filter(padjust <= 0.05) %>%
dplyr::select(Description, everything())
common_genes <- ATG7_res_ATG7_only_signif %>%
dplyr::count(Description) %>%
filter(n > 1) %>%
arrange(desc(n))
common_genes <-
ATG7_res_ATG7_only %>%
filter(Description %in% common_genes$Description,
sampleID %in% ATG7_patient_metadata[["samp_id_tidy"]]) %>%
arrange(Description) %>%
mutate(pval_zscore = str_c(zScore, ";", round(padjust, 2))) %>%
dplyr::select(`Sample ID` = sampleID, pval_zscore, gene_name = Description) %>%
arrange(gene_name) %>%
spread(key = gene_name, value = pval_zscore) %>%
left_join(ATG7_patient_metadata %>%
dplyr::select(samp_id_tidy, family),
by = c("Sample ID" = "samp_id_tidy")) %>%
dplyr::select(Family = family, everything()) %>%
arrange(Family)
common_genes %>% data.table::data.table()write_delim(common_genes, here::here("results/tables_figures/OUTRIDER_common_outliers.csv"), delim = ",")3.2.4 Exploratory analysis
3.2.4.1 Common outliers to M1856-17, M1111-19, and M1716-19
ATG7_res_ATG7_only %>%
filter(sampleID %in% c("M1856.17", "M1111.19", "M1716.19")) %>%
group_by(sampleID) %>%
mutate(rank = rank(padjust)) %>%
group_by(geneID) %>%
mutate(mean_rank = mean(rank)) %>%
ungroup() %>%
arrange(mean_rank) %>%
head(200) %>%
DT::datatable()# METTL18
plot_fpkm("ENSG00000171806", gene_count_ods, ATG7_patient_metadata)
# TMEM246
plot_fpkm("ENSG00000165152", gene_count_ods, ATG7_patient_metadata)
3.2.4.2 Negative controls
There are no outlier genes detected in either of the control samples, which is reassuring.
ATG7_res_all %>%
filter(str_detect(sampleID, "control")) %>%
mutate(padjust = p.adjust(pValue, method = "fdr")) %>%
filter(padjust <= 0.05)As a double check, how many of the ATG7 outliers are common to any of the remaining patient samples.
signif_genes_remaining_samp <- ATG7_res_all %>%
filter(!(sampleID %in% ATG7_patient_metadata[["samp_id_tidy"]])) %>%
mutate(padjust = p.adjust(pValue, method = "fdr")) %>%
filter(padjust <= 0.05) %>%
.[["Description"]]
mean(unique(ATG7_res_ATG7_only_signif$Description) %in% unique(signif_genes_remaining_samp))## [1] 0.20833333.2.4.3 Gene-set enrichment
Here, we use a ranked gprofiler test with all 264 outlier genes as input. The background genes used are all genes that have been tested through OUTRIDER.
hits_gprofiler <-
gprofiler2::gost(query = unique(ATG7_res_ATG7_only_signif %>%
arrange(padjust) %>%
.[["geneID"]]),
organism = "hsapiens",
ordered_query = T,
significant = T,
user_threshold = 0.05,
correction_method = c("g_SCS"),
domain_scope = c("custom"),
custom_bg = unique(ATG7_res_all$geneID),
sources = c("GO", "KEGG", "REAC"),
evcodes = F)
hits_gprofiler$result %>%
data.table::data.table()gprofiler2::gostplot(hits_gprofiler)3.2.4.4 Known autophagy hits
How many known autophagy genes that are part of the KEGG pathway are found amongst the ATG7 outliers? RAB39B comes up as upregulated in one patient (M0921.18). Though this is different from the RAB9 suggested by Jack.
query_kegg_pathway <- function(pathway_name, organism = "hsa"){
# get all pathways associated with humans
all_pathways <- keggList("pathway", organism)
# filter to only selected pathway and get information about it
selected_pathway <- all_pathways[str_detect(all_pathways, pathway_name)]
print(str_c("Getting info for the pathways: ", unname(selected_pathway) %>% str_c(collapse = ", ")))
selected_pathway_info <- keggGet(names(selected_pathway))
return(selected_pathway_info)
}
get_pathway_gene_info <- function(kegg_pathway_info, organism = "hsa"){
pathway_gene_info_all <- tibble()
for(i in 1:length(kegg_pathway_info)){
odds <- (1:length(kegg_pathway_info[[i]]$GENE))[(1:length(kegg_pathway_info[[i]]$GENE) %% 2) != 0]
evens <- (1:length(kegg_pathway_info[[i]]$GENE))[(1:length(kegg_pathway_info[[i]]$GENE) %% 2) == 0]
pathway_gene_info <-
tibble(pathway_kegg_code = kegg_pathway_info[[i]]$ENTRY,
pathway_name = kegg_pathway_info[[i]]$NAME,
organism = organism,
kegg_gene_entry = kegg_pathway_info[[i]]$GENE[odds],
gene = kegg_pathway_info[[i]]$GENE[evens]) %>%
separate(col = "gene", into = c("gene_name", "gene_desc"), sep = "; ")
pathway_gene_info_all <- pathway_gene_info_all %>% bind_rows(pathway_gene_info)
}
# for(i in 1:nrow(pathway_gene_info)){
#
# gene_info <- keggGet(str_c(organism, ":", pathway_gene_info$kegg_gene_entry)[i])
# pathway_gene_info$gene_identifiers[i] <- gene_info[[1]]$DBLINKS %>% str_c(collapse = ";")
#
# }
return(pathway_gene_info_all)
}
autophagy_kegg_info <- query_kegg_pathway(pathway_name = "Autophagy", organism = "hsa")## [1] "Getting info for the pathways: Autophagy - other - Homo sapiens (human), Autophagy - animal - Homo sapiens (human)"autophagy_gene_info <- get_pathway_gene_info(kegg_pathway_info = autophagy_kegg_info, organism = "hsa")
ATG7_res_ATG7_only_signif[ATG7_res_ATG7_only_signif$Description %in% unique(autophagy_gene_info$gene_name),]plot_fpkm(gene_id = "ENSG00000155961", gene_count_ods, ATG7_patient_metadata)
NFE2L2 may be up-regulated in ATG7 patients. NFE2L2 is not an expression outlier and just misses the p-value cut-off for DESeq2 (p-value: 0.07). However, when plotting the expression of NFE2L2 it looks like it may be up-regulated in ATG7 patients.
ATG7_res_ATG7_only %>% filter(Description == "NFE2L2")plot_fpkm(gene_id = "ENSG00000116044", gene_count_ods, ATG7_patient_metadata)
3.3 Mitosis and cell-cycle-related pathways are commonly disrupted across the ATG7 patients
DESeq2 detected 358 genes that were significantly differentially expressed with a log-fold change of >1.5 in the ATG7 patients (see full results in ATG7_deseq2_signif_logfold_ab_1.5.csv). The vast majority of these (210 genes) were down-regulated. When running gene-set enrichment analysis, we found that these 358 genes were generally enriched in mitotic and cell cycle related pathways (see full results in ATG7_deseq2_signif_logfold_ab_1.5_gprofiler.csv). However, four Golgi-related pathways were among those enriched amongst down-regulated genes (COPI-dependent Golgi-to-ER retrograde traffic, Golgi Cisternae Pericentriolar Stack Reorganization, Intra-Golgi and retrograde Golgi-to-ER traffic, Golgi-to-ER retrograde transport). We did not find that the differentially expressed genes were enriched for autophagy related pathways.
load(here::here("results/aberrant_expression/DESeq2/ATG7_deseq2_res.rda"))
ATG7_deseq2_res <- ATG7_deseq2_res %>%
filter(padj <= 0.05,
abs(log2FoldChange) >= 1.5) %>%
mutate(log2FoldChange = -log2FoldChange)
# print all genes
ATG7_deseq2_res %>%
arrange(log2FoldChange) %>%
DT::datatable()hits_gprofiler <-
gprofiler2::gost(query = ATG7_deseq2_res %>%
arrange(pvalue) %>%
.[["gene_id"]],
organism = "hsapiens",
ordered_query = TRUE,
significant = T,
user_threshold = 0.05,
correction_method = c("g_SCS"),
domain_scope = c("custom"),
custom_bg = unique(ATG7_res_all$geneID),
sources = c("GO", "KEGG", "REAC"),
evcodes = F)
# top 5 significant gprofiler hits
hits_gprofiler$result %>%
as_tibble() %>%
.[["term_name"]] %>% head()## [1] "cell cycle" "cell cycle process"
## [3] "mitotic cell cycle process" "cell division"
## [5] "mitotic cell cycle" "chromosome segregation"# golgi hits
hits_gprofiler$result %>%
as_tibble() %>%
filter(str_detect(term_name, "Golgi")) %>%
.[["term_name"]]## [1] "retrograde vesicle-mediated transport, Golgi to endoplasmic reticulum"
## [2] "COPI-dependent Golgi-to-ER retrograde traffic"
## [3] "Golgi-to-ER retrograde transport"
## [4] "Golgi Cisternae Pericentriolar Stack Reorganization"
## [5] "Intra-Golgi and retrograde Golgi-to-ER traffic"hits_gprofiler$result %>%
dplyr::select(-parents) %>%
write_delim(here::here("results/tmp/ATG7_deseq2_signif_logfold_ab_1.5_gprofiler.csv"), delim = ",")
gprofiler2::gostplot(hits_gprofiler)Although we found that mitotic-related pathways were commonly disrupted across ATG7 patients, it was not clear whether this was variable across the different ATG7 patients or families. In order to investigate this, we selected the gene PTPRD which was amongst the most down-regulated genes with the greatest negative log-fold change, as representative of the commonly disrupted mitotic pathways. The expression of PTPRD was lowest in family 1 and family 3. Whilst family 2/4 had relatively unperturbed PTPRD expression.
plot_fpkm(gene_id = "ENSG00000153707", gene_count_ods, ATG7_patient_metadata) +
scale_y_log10()
3.3.1 Exploratory analysis
3.3.1.1 Expression of other dysregulated genes from the mitotic pathway.
# NGFR
plot_fpkm(gene_id = "ENSG00000064300", gene_count_ods, ATG7_patient_metadata) 
# RARB
plot_fpkm(gene_id = "ENSG00000077092", gene_count_ods, ATG7_patient_metadata) 
# RSPO4
plot_fpkm(gene_id = "ENSG00000101282", gene_count_ods, ATG7_patient_metadata) 
# SOX11
plot_fpkm(gene_id = "ENSG00000176887", gene_count_ods, ATG7_patient_metadata) 
# TRIM55
plot_fpkm(gene_id = "ENSG00000147573", gene_count_ods, ATG7_patient_metadata) 
3.3.1.2 Known autophagy hits
ATG7_deseq2_res[ATG7_deseq2_res$gene_symbol %in% unique(autophagy_gene_info$gene_name),]# DEPTOR
plot_fpkm(gene_id = "ENSG00000155792", gene_count_ods, ATG7_patient_metadata) 
# DAPK1
plot_fpkm(gene_id = "ENSG00000196730", gene_count_ods, ATG7_patient_metadata) 
3.3.1.3 RAB9
ATG7_deseq2_res %>%
filter(gene_symbol == "RAB9A")ATG7_res_ATG7_only %>%
filter(Description == "RAB9A")plot_fpkm(gene_id = "ENSG00000123595", gene_count_ods, ATG7_patient_metadata)
3.3.1.4 Negative control
- The previous set of results were purely tagging the effect of batch variation. Here, I have also reun DESeq2 on the 2 COL6A samples as a negative control.
- The results, reassuringly, look to be completely different.
load(here::here("results/aberrant_expression/DESeq2/COL6A_deseq2_res.rda"))
hits_gprofiler <-
gprofiler2::gost(query = COL6A_deseq2_res %>%
filter(padj <= 0.05, abs(log2FoldChange) >= 1.5) %>%
arrange(pvalue) %>%
.[["gene_id"]],
organism = "hsapiens",
ordered_query = T,
significant = T,
user_threshold = 0.05,
correction_method = c("g_SCS"),
domain_scope = c("custom"),
custom_bg = unique(ATG7_res_all$geneID),
sources = c("GO", "KEGG", "REAC"),
evcodes = F)
hits_gprofiler$result %>%
dplyr::select(-parents) %>%
write_delim(here::here("results/tmp/COL6A_deseq2_signif_logfold_ab_1.5_gprofiler_all.csv"), delim = ",")
hits_gprofiler$result %>%
data.table::data.table()gprofiler2::gostplot(hits_gprofiler)4 Reproducibility
## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
## setting value
## version R version 3.6.1 (2019-07-05)
## os Ubuntu 16.04.6 LTS
## system x86_64, linux-gnu
## ui X11
## language (EN)
## collate en_GB.UTF-8
## ctype en_GB.UTF-8
## tz Europe/London
## date 2020-10-15
##
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
## package * version date lib source
## abind 1.4-5 2016-07-21 [2] CRAN (R 3.6.1)
## annotate 1.64.0 2019-10-29 [1] Bioconductor
## AnnotationDbi * 1.48.0 2019-10-29 [1] Bioconductor
## askpass 1.1 2019-01-13 [2] CRAN (R 3.6.1)
## assertthat 0.2.1 2019-03-21 [2] CRAN (R 3.6.1)
## backports 1.1.10 2020-09-15 [1] CRAN (R 3.6.1)
## base64enc 0.1-3 2015-07-28 [2] CRAN (R 3.6.1)
## BBmisc 1.11 2017-03-10 [1] CRAN (R 3.6.1)
## Biobase * 2.46.0 2019-10-29 [1] Bioconductor
## BiocFileCache 1.10.2 2019-11-08 [1] Bioconductor
## BiocGenerics * 0.32.0 2019-10-29 [1] Bioconductor
## BiocParallel * 1.20.1 2019-12-21 [1] Bioconductor
## biomaRt 2.42.1 2020-03-26 [1] Bioconductor
## Biostrings 2.54.0 2019-10-29 [1] Bioconductor
## bit 4.0.4 2020-08-04 [1] CRAN (R 3.6.1)
## bit64 4.0.5 2020-08-30 [1] CRAN (R 3.6.1)
## bitops 1.0-6 2013-08-17 [2] CRAN (R 3.6.1)
## blob 1.2.1 2020-01-20 [1] CRAN (R 3.6.1)
## broom 0.5.6 2020-04-20 [1] CRAN (R 3.6.1)
## car 3.0-8 2020-05-21 [1] CRAN (R 3.6.1)
## carData 3.0-4 2020-05-22 [1] CRAN (R 3.6.1)
## caTools 1.18.0 2020-01-17 [1] CRAN (R 3.6.1)
## cellranger 1.1.0 2016-07-27 [2] CRAN (R 3.6.1)
## checkmate 2.0.0 2020-02-06 [1] CRAN (R 3.6.1)
## cli 2.1.0 2020-10-12 [1] CRAN (R 3.6.1)
## cluster 2.1.0 2019-06-19 [4] CRAN (R 3.6.1)
## codetools 0.2-16 2018-12-24 [4] CRAN (R 3.6.1)
## colorspace 1.4-1 2019-03-18 [2] CRAN (R 3.6.1)
## crayon 1.3.4 2017-09-16 [2] CRAN (R 3.6.1)
## crosstalk 1.1.0.1 2020-03-13 [1] CRAN (R 3.6.1)
## curl 4.3 2019-12-02 [1] CRAN (R 3.6.1)
## data.table * 1.13.0 2020-07-24 [1] CRAN (R 3.6.1)
## DBI 1.1.0 2019-12-15 [2] CRAN (R 3.6.1)
## dbplyr 1.4.4 2020-05-27 [1] CRAN (R 3.6.1)
## DelayedArray * 0.12.3 2020-04-09 [1] Bioconductor
## dendextend 1.14.0 2020-08-26 [1] CRAN (R 3.6.1)
## DESeq2 1.26.0 2019-10-29 [1] Bioconductor
## digest 0.6.25 2020-02-23 [1] CRAN (R 3.6.1)
## dplyr * 1.0.2 2020-08-18 [1] CRAN (R 3.6.1)
## DT 0.14 2020-06-24 [1] CRAN (R 3.6.1)
## ellipsis 0.3.1 2020-05-15 [1] CRAN (R 3.6.1)
## evaluate 0.14 2019-05-28 [2] CRAN (R 3.6.1)
## fansi 0.4.1 2020-01-08 [1] CRAN (R 3.6.1)
## farver 2.0.3 2020-01-16 [1] CRAN (R 3.6.1)
## fastmap 1.0.1 2019-10-08 [1] CRAN (R 3.6.1)
## forcats * 0.5.0 2020-03-01 [2] CRAN (R 3.6.1)
## foreach 1.5.0 2020-03-30 [1] CRAN (R 3.6.1)
## foreign 0.8-72 2019-08-02 [4] CRAN (R 3.6.1)
## Formula 1.2-3 2018-05-03 [2] CRAN (R 3.6.1)
## fs 1.5.0 2020-07-31 [1] CRAN (R 3.6.1)
## gclus 1.3.2 2019-01-07 [1] CRAN (R 3.6.1)
## gdata 2.18.0 2017-06-06 [1] CRAN (R 3.6.1)
## genefilter 1.68.0 2019-10-29 [1] Bioconductor
## geneplotter 1.64.0 2019-10-29 [1] Bioconductor
## generics 0.0.2 2018-11-29 [1] CRAN (R 3.6.1)
## GenomeInfoDb * 1.22.1 2020-03-27 [1] Bioconductor
## GenomeInfoDbData 1.2.2 2020-08-13 [1] Bioconductor
## GenomicAlignments 1.22.1 2019-11-12 [1] Bioconductor
## GenomicFeatures * 1.38.2 2020-02-15 [1] Bioconductor
## GenomicRanges * 1.38.0 2019-10-29 [1] Bioconductor
## ggplot2 * 3.3.2 2020-06-19 [1] CRAN (R 3.6.1)
## ggpubr 0.4.0 2020-06-27 [1] CRAN (R 3.6.1)
## ggsci 2.9 2018-05-14 [1] CRAN (R 3.6.1)
## ggsignif 0.6.0 2019-08-08 [1] CRAN (R 3.6.1)
## glue 1.4.2 2020-08-27 [1] CRAN (R 3.6.1)
## gplots 3.0.4 2020-07-05 [1] CRAN (R 3.6.1)
## gprofiler2 0.1.9 2020-04-23 [1] CRAN (R 3.6.1)
## gridExtra 2.3 2017-09-09 [2] CRAN (R 3.6.1)
## gtable 0.3.0 2019-03-25 [2] CRAN (R 3.6.1)
## gtools 3.8.2 2020-03-31 [1] CRAN (R 3.6.1)
## haven 2.3.1 2020-06-01 [1] CRAN (R 3.6.1)
## heatmaply 1.1.1 2020-08-25 [1] CRAN (R 3.6.1)
## here 0.1 2017-05-28 [1] CRAN (R 3.6.1)
## Hmisc 4.4-1 2020-08-10 [1] CRAN (R 3.6.1)
## hms 0.5.3 2020-01-08 [1] CRAN (R 3.6.1)
## htmlTable 2.1.0 2020-09-16 [1] CRAN (R 3.6.1)
## htmltools 0.5.0 2020-06-16 [1] CRAN (R 3.6.1)
## htmlwidgets 1.5.1 2019-10-08 [1] CRAN (R 3.6.1)
## httpuv 1.5.2 2019-09-11 [2] CRAN (R 3.6.1)
## httr 1.4.2 2020-07-20 [1] CRAN (R 3.6.1)
## IRanges * 2.20.2 2020-01-13 [1] Bioconductor
## iterators 1.0.12 2019-07-26 [2] CRAN (R 3.6.1)
## jpeg 0.1-8.1 2019-10-24 [1] CRAN (R 3.6.1)
## jsonlite 1.7.1 2020-09-07 [1] CRAN (R 3.6.1)
## KEGGREST * 1.24.1 2019-10-08 [1] Bioconductor
## KernSmooth 2.23-17 2020-04-26 [1] CRAN (R 3.6.1)
## knitr 1.29 2020-06-23 [1] CRAN (R 3.6.1)
## labeling 0.3 2014-08-23 [2] CRAN (R 3.6.1)
## later 1.1.0.1 2020-06-05 [1] CRAN (R 3.6.1)
## lattice 0.20-38 2018-11-04 [4] CRAN (R 3.6.1)
## latticeExtra 0.6-29 2019-12-19 [2] CRAN (R 3.6.1)
## lazyeval 0.2.2 2019-03-15 [2] CRAN (R 3.6.1)
## lifecycle 0.2.0 2020-03-06 [1] CRAN (R 3.6.1)
## locfit 1.5-9.4 2020-03-25 [1] CRAN (R 3.6.1)
## lubridate 1.7.9 2020-06-08 [1] CRAN (R 3.6.1)
## magrittr 1.5 2014-11-22 [2] CRAN (R 3.6.1)
## MASS 7.3-51.6 2020-04-26 [1] CRAN (R 3.6.1)
## Matrix 1.2-18 2019-11-27 [1] CRAN (R 3.6.1)
## matrixStats * 0.57.0 2020-09-25 [1] CRAN (R 3.6.1)
## memoise 1.1.0 2017-04-21 [2] CRAN (R 3.6.1)
## mime 0.9 2020-02-04 [1] CRAN (R 3.6.1)
## modelr 0.1.5 2019-08-08 [2] CRAN (R 3.6.1)
## munsell 0.5.0 2018-06-12 [2] CRAN (R 3.6.1)
## nlme 3.1-148 2020-05-24 [1] CRAN (R 3.6.1)
## nnet 7.3-12 2016-02-02 [4] CRAN (R 3.6.1)
## openssl 1.4.3 2020-09-18 [1] CRAN (R 3.6.1)
## openxlsx 4.1.5 2020-05-06 [1] CRAN (R 3.6.1)
## OUTRIDER * 1.5.1 2020-09-17 [1] Github (gagneurlab/OUTRIDER@7fbdaaf)
## pcaMethods 1.78.0 2019-10-29 [1] Bioconductor
## pheatmap 1.0.12 2019-01-04 [1] CRAN (R 3.6.1)
## pillar 1.4.6 2020-07-10 [1] CRAN (R 3.6.1)
## pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 3.6.1)
## plotly 4.9.2.1 2020-04-04 [1] CRAN (R 3.6.1)
## plyr 1.8.6 2020-03-03 [2] CRAN (R 3.6.1)
## png 0.1-7 2013-12-03 [1] CRAN (R 3.6.1)
## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 3.6.1)
## progress 1.2.2 2019-05-16 [2] CRAN (R 3.6.1)
## promises 1.1.1 2020-06-09 [1] CRAN (R 3.6.1)
## PRROC 1.3.1 2018-06-19 [1] CRAN (R 3.6.1)
## purrr * 0.3.4 2020-04-17 [1] CRAN (R 3.6.1)
## R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.1)
## rappdirs 0.3.1 2016-03-28 [1] CRAN (R 3.6.1)
## RColorBrewer 1.1-2 2014-12-07 [2] CRAN (R 3.6.1)
## Rcpp 1.0.5 2020-07-06 [1] CRAN (R 3.6.1)
## RCurl 1.98-1.2 2020-04-18 [1] CRAN (R 3.6.1)
## readr * 1.4.0 2020-10-05 [1] CRAN (R 3.6.1)
## readxl 1.3.1 2019-03-13 [2] CRAN (R 3.6.1)
## registry 0.5-1 2019-03-05 [1] CRAN (R 3.6.1)
## reprex 0.3.0 2019-05-16 [2] CRAN (R 3.6.1)
## reshape2 1.4.4 2020-04-09 [2] CRAN (R 3.6.1)
## rio 0.5.16 2018-11-26 [1] CRAN (R 3.6.1)
## rlang 0.4.8 2020-10-08 [1] CRAN (R 3.6.1)
## rmarkdown 2.3 2020-06-18 [1] CRAN (R 3.6.1)
## rpart 4.1-15 2019-04-12 [4] CRAN (R 3.6.1)
## rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.6.1)
## Rsamtools 2.2.3 2020-02-23 [1] Bioconductor
## RSQLite 2.2.1 2020-09-30 [1] CRAN (R 3.6.1)
## rstatix 0.6.0 2020-06-18 [1] CRAN (R 3.6.1)
## rstudioapi 0.11 2020-02-07 [2] CRAN (R 3.6.1)
## rtracklayer 1.46.0 2019-10-29 [1] Bioconductor
## rvest 0.3.5 2019-11-08 [1] CRAN (R 3.6.1)
## S4Vectors * 0.24.4 2020-04-09 [1] Bioconductor
## scales 1.1.1 2020-05-11 [1] CRAN (R 3.6.1)
## seriation 1.2-8 2019-08-27 [1] CRAN (R 3.6.1)
## sessioninfo * 1.1.1 2018-11-05 [2] CRAN (R 3.6.1)
## shiny 1.5.0 2020-06-23 [1] CRAN (R 3.6.1)
## stringi 1.5.3 2020-09-09 [1] CRAN (R 3.6.1)
## stringr * 1.4.0 2019-02-10 [2] CRAN (R 3.6.1)
## SummarizedExperiment * 1.16.1 2019-12-19 [1] Bioconductor
## survival 3.2-3 2020-06-13 [1] CRAN (R 3.6.1)
## tibble * 3.0.4 2020-10-12 [1] CRAN (R 3.6.1)
## tidyr * 1.1.2 2020-08-27 [1] CRAN (R 3.6.1)
## tidyselect 1.1.0 2020-05-11 [1] CRAN (R 3.6.1)
## tidyverse * 1.3.0 2019-11-21 [1] CRAN (R 3.6.1)
## TSP 1.1-10 2020-04-17 [1] CRAN (R 3.6.1)
## vctrs 0.3.4 2020-08-29 [1] CRAN (R 3.6.1)
## viridis 0.5.1 2018-03-29 [2] CRAN (R 3.6.1)
## viridisLite 0.3.0 2018-02-01 [2] CRAN (R 3.6.1)
## webshot 0.5.2 2019-11-22 [1] CRAN (R 3.6.1)
## withr 2.3.0 2020-09-22 [1] CRAN (R 3.6.1)
## xfun 0.18 2020-09-29 [1] CRAN (R 3.6.1)
## XML 3.99-0.3 2020-01-20 [2] CRAN (R 3.6.1)
## xml2 1.3.2 2020-04-23 [1] CRAN (R 3.6.1)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 3.6.1)
## XVector 0.26.0 2019-10-29 [1] Bioconductor
## yaml 2.2.1 2020-02-01 [1] CRAN (R 3.6.1)
## zip 2.0.4 2019-09-01 [1] CRAN (R 3.6.1)
## zlibbioc 1.32.0 2019-10-29 [1] Bioconductor
##
## [1] /home/dzhang/R/x86_64-pc-linux-gnu-library/3.6
## [2] /usr/local/lib/R/site-library
## [3] /usr/lib/R/site-library
## [4] /usr/lib/R/library